home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / LINUX / MATH_EMU.ZIP / MATH_EMU / FPU_SQRT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1979-12-31  |  14.7 KB  |  392 lines
  1. /*        $NetBSD: fpu_sqrt.c,v 1.2 1994/11/20 20:52:46 deraadt Exp $ */
  2.  
  3. /*
  4.  * Copyright (c) 1992, 1993
  5.  *        The Regents of the University of California.  All rights reserved.
  6.  *
  7.  * This software was developed by the Computer Systems Engineering group
  8.  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
  9.  * contributed to Berkeley.
  10.  *
  11.  * All advertising materials mentioning features or use of this software
  12.  * must display the following acknowledgement:
  13.  *        This product includes software developed by the University of
  14.  *        California, Lawrence Berkeley Laboratory.
  15.  *
  16.  * Redistribution and use in source and binary forms, with or without
  17.  * modification, are permitted provided that the following conditions
  18.  * are met:
  19.  * 1. Redistributions of source code must retain the above copyright
  20.  *    notice, this list of conditions and the following disclaimer.
  21.  * 2. Redistributions in binary form must reproduce the above copyright
  22.  *    notice, this list of conditions and the following disclaimer in the
  23.  *    documentation and/or other materials provided with the distribution.
  24.  * 3. All advertising materials mentioning features or use of this software
  25.  *    must display the following acknowledgement:
  26.  *        This product includes software developed by the University of
  27.  *        California, Berkeley and its contributors.
  28.  * 4. Neither the name of the University nor the names of its contributors
  29.  *    may be used to endorse or promote products derived from this software
  30.  *    without specific prior written permission.
  31.  *
  32.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  33.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  34.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  35.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  36.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  37.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  38.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  39.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  40.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  41.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  42.  * SUCH DAMAGE.
  43.  *
  44.  *        @(#)fpu_sqrt.c      8.1 (Berkeley) 6/11/93
  45.  */
  46.  
  47. /*
  48.  * Perform an FPU square root (return sqrt(x)).
  49.  */
  50.  
  51. #include "types.h"
  52.  
  53. #include "reg.h"
  54. #include "fpu_arit.h"
  55. #include "fpu_emul.h"
  56.  
  57. /*
  58.  * Our task is to calculate the square root of a floating point number x0.
  59.  * This number x normally has the form:
  60.  *
  61.  *                      exp
  62.  *        x = mant * 2                  (where 1 <= mant < 2 and exp is an integer)
  63.  *
  64.  * This can be left as it stands, or the mantissa can be doubled and the
  65.  * exponent decremented:
  66.  *
  67.  *                              exp-1
  68.  *        x = (2 * mant) * 2  (where 2 <= 2 * mant < 4)
  69.  *
  70.  * If the exponent `exp' is even, the square root of the number is best
  71.  * handled using the first form, and is by definition equal to:
  72.  *
  73.  *                                      exp/2
  74.  *        sqrt(x) = sqrt(mant) * 2
  75.  *
  76.  * If exp is odd, on the other hand, it is convenient to use the second
  77.  * form, giving:
  78.  *
  79.  *                                          (exp-1)/2
  80.  *        sqrt(x) = sqrt(2 * mant) * 2
  81.  *
  82.  * In the first case, we have
  83.  *
  84.  *        1 <= mant < 2
  85.  *
  86.  * and therefore
  87.  *
  88.  *        sqrt(1) <= sqrt(mant) < sqrt(2)
  89.  *
  90.  * while in the second case we have
  91.  *
  92.  *        2 <= 2*mant < 4
  93.  *
  94.  * and therefore
  95.  *
  96.  *        sqrt(2) <= sqrt(2*mant) < sqrt(4)
  97.  *
  98.  * so that in any case, we are sure that
  99.  *
  100.  *        sqrt(1) <= sqrt(n * mant) < sqrt(4),    n = 1 or 2
  101.  *
  102.  * or
  103.  *
  104.  *        1 <= sqrt(n * mant) < 2,                n = 1 or 2.
  105.  *
  106.  * This root is therefore a properly formed mantissa for a floating
  107.  * point number.  The exponent of sqrt(x) is either exp/2 or (exp-1)/2
  108.  * as above.  This leaves us with the problem of finding the square root
  109.  * of a fixed-point number in the range [1..4).
  110.  *
  111.  * Though it may not be instantly obvious, the following square root
  112.  * algorithm works for any integer x of an even number of bits, provided
  113.  * that no overflows occur:
  114.  *
  115.  *        let q = 0
  116.  *        for k = NBITS-1 to 0 step -1 do -- for each digit in the answer...
  117.  *                  x *= 2                        -- multiply by radix, for next digit
  118.  *                  if x >= 2q + 2^k then         -- if adding 2^k does not
  119.  *                            x -= 2q + 2^k       -- exceed the correct root,
  120.  *                            q += 2^k  -- add 2^k and adjust x
  121.  *                  fi
  122.  *        done
  123.  *        sqrt = q / 2^(NBITS/2)                  -- (and any remainder is in x)
  124.  *
  125.  * If NBITS is odd (so that k is initially even), we can just add another
  126.  * zero bit at the top of x.  Doing so means that q is not going to acquire
  127.  * a 1 bit in the first trip around the loop (since x0 < 2^NBITS).  If the
  128.  * final value in x is not needed, or can be off by a factor of 2, this is
  129.  * equivalant to moving the `x *= 2' step to the bottom of the loop:
  130.  *
  131.  *        for k = NBITS-1 to 0 step -1 do if ... fi; x *= 2; done
  132.  *
  133.  * and the result q will then be sqrt(x0) * 2^floor(NBITS / 2).
  134.  * (Since the algorithm is destructive on x, we will call x's initial
  135.  * value, for which q is some power of two times its square root, x0.)
  136.  *
  137.  * If we insert a loop invariant y = 2q, we can then rewrite this using
  138.  * C notation as:
  139.  *
  140.  *        q = y = 0; x = x0;
  141.  *        for (k = NBITS; --k >= 0;) {
  142.  * #if (NBITS is even)
  143.  *                  x *= 2;
  144.  * #endif
  145.  *                  t = y + (1 << k);
  146.  *                  if (x >= t) {
  147.  *                            x -= t;
  148.  *                            q += 1 << k;
  149.  *                            y += 1 << (k + 1);
  150.  *                  }
  151.  * #if (NBITS is odd)
  152.  *                  x *= 2;
  153.  * #endif
  154.  *        }
  155.  *
  156.  * If x0 is fixed point, rather than an integer, we can simply alter the
  157.  * scale factor between q and sqrt(x0).  As it happens, we can easily arrange
  158.  * for the scale factor to be 2**0 or 1, so that sqrt(x0) == q.
  159.  *
  160.  * In our case, however, x0 (and therefore x, y, q, and t) are multiword
  161.  * integers, which adds some complication.  But note that q is built one
  162.  * bit at a time, from the top down, and is not used itself in the loop
  163.  * (we use 2q as held in y instead).  This means we can build our answer
  164.  * in an integer, one word at a time, which saves a bit of work.  Also,
  165.  * since 1 << k is always a `new' bit in q, 1 << k and 1 << (k+1) are
  166.  * `new' bits in y and we can set them with an `or' operation rather than
  167.  * a full-blown multiword add.
  168.  *
  169.  * We are almost done, except for one snag.  We must prove that none of our
  170.  * intermediate calculations can overflow.  We know that x0 is in [1..4)
  171.  * and therefore the square root in q will be in [1..2), but what about x,
  172.  * y, and t?
  173.  *
  174.  * We know that y = 2q at the beginning of each loop.  (The relation only
  175.  * fails temporarily while y and q are being updated.)  Since q < 2, y < 4.
  176.  * The sum in t can, in our case, be as much as y+(1<<1) = y+2 < 6, and.
  177.  * Furthermore, we can prove with a bit of work that x never exceeds y by
  178.  * more than 2, so that even after doubling, 0 <= x < 8.  (This is left as
  179.  * an exercise to the reader, mostly because I have become tired of working
  180.  * on this comment.)
  181.  *
  182.  * If our floating point mantissas (which are of the form 1.frac) occupy
  183.  * B+1 bits, our largest intermediary needs at most B+3 bits, or two extra.
  184.  * In fact, we want even one more bit (for a carry, to avoid compares), or
  185.  * three extra.  There is a comment in fpu_emu.h reminding maintainers of
  186.  * this, so we have some justification in assuming it.
  187.  */
  188. struct fpn *
  189. fpu_sqrt(fe)
  190.           struct fpemu *fe;
  191. {
  192.           register struct fpn *x = &fe->fe_f2;
  193.           register u_int bit, q, tt;
  194.           register u_int x0, x1, x2, x3;
  195.           register u_int y0, y1, y2, y3;
  196.           register u_int d0, d1, d2, d3;
  197.           register int e;
  198.           FPU_DECL_CARRY
  199.  
  200.           /*
  201.            * Take care of special cases first.  In order:
  202.            *
  203.            *        sqrt(NaN) = NaN
  204.            *        sqrt(+0) = +0
  205.            *        sqrt(-0) = -0
  206.            *        sqrt(x < 0) = NaN   (including sqrt(-Inf))
  207.            *        sqrt(+Inf) = +Inf
  208.            *
  209.            * Then all that remains are numbers with mantissas in [1..2).
  210.            */
  211.           if (ISNAN(x) || ISZERO(x))
  212.                     return (x);
  213.           if (x->fp_sign)
  214.                     return (fpu_newnan(fe));
  215.           if (ISINF(x))
  216.                     return (x);
  217.  
  218.           /*
  219.            * Calculate result exponent.  As noted above, this may involve
  220.            * doubling the mantissa.  We will also need to double x each
  221.            * time around the loop, so we define a macro for this here, and
  222.            * we break out the multiword mantissa.
  223.            */
  224. #ifdef FPU_SHL1_BY_ADD
  225. #define   DOUBLE_X { FPU_ADDS(x3, x3, x3); FPU_ADDCS(x2, x2, x2); FPU_ADDCS(x1, x1, x1); FPU_ADDC(x0, x0, x0); }
  226. #else
  227. #define   DOUBLE_X { x0 = (x0 << 1) | (x1 >> 31); x1 = (x1 << 1) | (x2 >> 31); x2 = (x2 << 1) | (x3 >> 31); x3 <<= 1; }
  228. #endif
  229. #if (FP_NMANT & 1) != 0
  230. # define ODD_DOUBLE DOUBLE_X
  231. # define EVEN_DOUBLE          /* nothing */
  232. #else
  233. # define ODD_DOUBLE /* nothing */
  234. # define EVEN_DOUBLE          DOUBLE_X
  235. #endif
  236.           x0 = x->fp_mant[0];
  237.           x1 = x->fp_mant[1];
  238.           x2 = x->fp_mant[2];
  239.           x3 = x->fp_mant[3];
  240.           e = x->fp_exp;
  241.           if (e & 1)                    /* exponent is odd; use sqrt(2mant) */
  242.                     DOUBLE_X;
  243.           /* THE FOLLOWING ASSUMES THAT RIGHT SHIFT DOES SIGN EXTENSION */
  244.           x->fp_exp = e >> 1; /* calculates (e&1 ? (e-1)/2 : e/2 */
  245.  
  246.           /*
  247.            * Now calculate the mantissa root.  Since x is now in [1..4),
  248.            * we know that the first trip around the loop will definitely
  249.            * set the top bit in q, so we can do that manually and start
  250.            * the loop at the next bit down instead.  We must be sure to
  251.            * double x correctly while doing the `known q=1.0'.
  252.            *
  253.            * We do this one mantissa-word at a time, as noted above, to
  254.            * save work.  To avoid `(1 << 31) << 1', we also do the top bit
  255.            * outside of each per-word loop.
  256.            *
  257.            * The calculation `t = y + bit' breaks down into `t0 = y0, ...,
  258.            * t3 = y3, t? |= bit' for the appropriate word.  Since the bit
  259.            * is always a `new' one, this means that three of the `t?'s are
  260.            * just the corresponding `y?'; we use `#define's here for this.
  261.            * The variable `tt' holds the actual `t?' variable.
  262.            */
  263.  
  264.           /* calculate q0 */
  265. #define   t0 tt
  266.           bit = FP_1;
  267.           EVEN_DOUBLE;
  268.           /* if (x >= (t0 = y0 | bit)) { */       /* always true */
  269.                     q = bit;
  270.                     x0 -= bit;
  271.                     y0 = bit << 1;
  272.           /* } */
  273.           ODD_DOUBLE;
  274.           while ((bit >>= 1) != 0) {    /* for remaining bits in q0 */
  275.                     EVEN_DOUBLE;
  276.                     t0 = y0 | bit;                /* t = y + bit */
  277.                     if (x0 >= t0) {               /* if x >= t then */
  278.                               x0 -= t0; /*        x -= t */
  279.                               q |= bit; /*        q += bit */
  280.                               y0 |= bit << 1;     /*        y += bit << 1 */
  281.                     }
  282.                     ODD_DOUBLE;
  283.           }
  284.           x->fp_mant[0] = q;
  285. #undef t0
  286.  
  287.           /* calculate q1.  note (y0&1)==0. */
  288. #define t0 y0
  289. #define t1 tt
  290.           q = 0;
  291.           y1 = 0;
  292.           bit = 1 << 31;
  293.           EVEN_DOUBLE;
  294.           t1 = bit;
  295.           FPU_SUBS(d1, x1, t1);
  296.           FPU_SUBC(d0, x0, t0);                   /* d = x - t */
  297.           if ((int)d0 >= 0) {           /* if d >= 0 (i.e., x >= t) then */
  298.                     x0 = d0, x1 = d1;   /*        x -= t */
  299.                     q = bit;            /*        q += bit */
  300.                     y0 |= 1;            /*        y += bit << 1 */
  301.           }
  302.           ODD_DOUBLE;
  303.           while ((bit >>= 1) != 0) {    /* for remaining bits in q1 */
  304.                     EVEN_DOUBLE;                  /* as before */
  305.                     t1 = y1 | bit;
  306.                     FPU_SUBS(d1, x1, t1);
  307.                     FPU_SUBC(d0, x0, t0);
  308.                     if ((int)d0 >= 0) {
  309.                               x0 = d0, x1 = d1;
  310.                               q |= bit;
  311.                               y1 |= bit << 1;
  312.                     }
  313.                     ODD_DOUBLE;
  314.           }
  315.           x->fp_mant[1] = q;
  316. #undef t1
  317.  
  318.           /* calculate q2.  note (y1&1)==0; y0 (aka t0) is fixed. */
  319. #define t1 y1
  320. #define t2 tt
  321.           q = 0;
  322.           y2 = 0;
  323.           bit = 1 << 31;
  324.           EVEN_DOUBLE;
  325.           t2 = bit;
  326.           FPU_SUBS(d2, x2, t2);
  327.           FPU_SUBCS(d1, x1, t1);
  328.           FPU_SUBC(d0, x0, t0);
  329.           if ((int)d0 >= 0) {
  330.                     x0 = d0, x1 = d1, x2 = d2;
  331.                     q |= bit;
  332.                     y1 |= 1;            /* now t1, y1 are set in concrete */
  333.           }
  334.           ODD_DOUBLE;
  335.           while ((bit >>= 1) != 0) {
  336.                     EVEN_DOUBLE;
  337.                     t2 = y2 | bit;
  338.                     FPU_SUBS(d2, x2, t2);
  339.                     FPU_SUBCS(d1, x1, t1);
  340.                     FPU_SUBC(d0, x0, t0);
  341.                     if ((int)d0 >= 0) {
  342.                               x0 = d0, x1 = d1, x2 = d2;
  343.                               q |= bit;
  344.                               y2 |= bit << 1;
  345.                     }
  346.                     ODD_DOUBLE;
  347.           }
  348.           x->fp_mant[2] = q;
  349. #undef t2
  350.  
  351.           /* calculate q3.  y0, t0, y1, t1 all fixed; y2, t2, almost done. */
  352. #define t2 y2
  353. #define t3 tt
  354.           q = 0;
  355.           y3 = 0;
  356.           bit = 1 << 31;
  357.           EVEN_DOUBLE;
  358.           t3 = bit;
  359.           FPU_SUBS(d3, x3, t3);
  360.           FPU_SUBCS(d2, x2, t2);
  361.           FPU_SUBCS(d1, x1, t1);
  362.           FPU_SUBC(d0, x0, t0);
  363.           ODD_DOUBLE;
  364.           if ((int)d0 >= 0) {
  365.                     x0 = d0, x1 = d1, x2 = d2;
  366.                     q |= bit;
  367.                     y2 |= 1;
  368.           }
  369.           while ((bit >>= 1) != 0) {
  370.                     EVEN_DOUBLE;
  371.                     t3 = y3 | bit;
  372.                     FPU_SUBS(d3, x3, t3);
  373.                     FPU_SUBCS(d2, x2, t2);
  374.                     FPU_SUBCS(d1, x1, t1);
  375.                     FPU_SUBC(d0, x0, t0);
  376.                     if ((int)d0 >= 0) {
  377.                               x0 = d0, x1 = d1, x2 = d2;
  378.                               q |= bit;
  379.                               y3 |= bit << 1;
  380.                     }
  381.                     ODD_DOUBLE;
  382.           }
  383.           x->fp_mant[3] = q;
  384.  
  385.           /*
  386.            * The result, which includes guard and round bits, is exact iff
  387.            * x is now zero; any nonzero bits in x represent sticky bits.
  388.            */
  389.           x->fp_sticky = x0 | x1 | x2 | x3;
  390.           return (x);
  391. }
  392.